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ABSTRACT 

We present the study of a set of N-body+SPH simulations of a Milky Way-like system 
produced by the radiative cooling of hot gas embedded in a dark matter halo. The 
galaxy and its gaseous halo evolve for 10 Gyr in isolation, which allows us to study how 
internal processes affect the evolution of the system. We show how the morphology, 
the kinematics and the evolution of the galaxy are affected by the input supernova 
feedback energy ^sn, and we compare its properties with those of the Milky Way. 
Different values of ^sn do not significantly affect the star formation history of the 
system, but the disc of cold gas gets thicker and more turbulent as feedback increases. 
Our main result is that, for the highest value of £’sn considered, the galaxy shows a 
prominent layer of extra-planar cold (log(T/K) < 4.3) gas extended up to a few kpc 
above the disc at column densities of 10^^ cm“^. The kinematics of this material is in 
agreement with that inferred for the Hi halos of our Galaxy and NGC 891, although 
its mass is lower. Also, the location, the kinematics and the typical column densities 
of the hot (5.3 < log(T/K) < 5.7) gas are in good agreement with those determined 
from the O VI absorption systems in the halo of the Milky Way and external galaxies. 
In contrast with the observations, however, gas at log(T/K) < 5.3 is lacking in the 
circumgalactic region of our systems. 

Key words: methods: numerical - Galaxy: evolution - Galaxy: halo - ISM: kine¬ 
matics and dynamics 


1 INTRODUCTION 


Numerical simulations constitute a fundamental tool to un¬ 
derstand the processes of galaxy formation and evolution. 
Simulations of isolated systems have proved to be very pow¬ 
erful to investigate the physics on galactic scales, and they 
have been widely used in the last decade to interpret the sec¬ 
ular evolution of discs (e.g. Debattista et al.||2006 Combes 


2006 


2008| , to study the role of stellar feedback (e.g. Stinson et al. 
Dalla Vecchia & Schaye 2008[ |2012| ), and to under¬ 


stand the effect of major and minor mergers on the mor¬ 
phology and dynamics of galaxies (e.g. 

2009| [Bois et al.|2011 ). 


Villalobos & Helmi 


Simulations can be used to address a long-standing 
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question in the theory of galaxy evolution, namely how disc 
galaxies are able to sustain star formation for a Hubble time 
without consuming their gas reservoir. In our Galaxy, for in¬ 
stance, the star formation rate (SFR) of the solar neighbor¬ 
hood has remained almost constant for the last ~ 10 Gyr 
(Twarog 1980| Aumer Binney||2QQ9| , implying that the 
gas consumed by star formation is continuously replenished. 
More generally, there is evidence that galaxies, at all times, 
must have accreted gas at a rate proportional to their SFR 
density ( Hopkins et al.|2d08 Fraternali Tomassetti|2012 ). 
However, there is little observational evidence for cold (Hi) 
gas accretion occurring onto galaxies at the required rate 
(e.g. Sancisi et al.||20Q8 ). 


ACDM cosmological simulations made a major break¬ 
through in this field by revealing the modes by which galax¬ 
ies may accrete their gas from the cosmic web. Low-mass 
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halos are expected to accrete material through the so-called 
‘cold-mode’, where filaments of gas at relatively low tem¬ 
perature (< 10^ K) are able to penetrate the halos down 
to the centre and feed directly the formation of the central 
galaxy (e.g. Keres et ah 2009). In massive halos, instead, 
accreting hlaments get shock-heated and settle into a rar¬ 
efied atmosphere of plasma at the virial temperature of the 
halo, often called coronae ( White fc Rees||1978 Crain et al. 


2010). Because of the hierarchical formation of structures. 


the high-redshift universe is expected to be dominated by 
the cold-mode accretion, which is progressively replaced by 
the ‘hot-mode’ as structures grow (Keres et al. p005l|DekeI| 


fc Birnboim||2006|. Massive disc galaxies at redshift z = 0 


including the Milky Way - are therefore expected to be 
surrounded by a large amount of rarefied hot gas - the so- 
called ‘coronae’ - which would constitute a vast reservoir of 
material to fuel star formation. 

Simulations also reveal that feedback from stars and 
active galactic nuclei (AGN) plays a fundamental role in 
the evolution of galaxies. In most cases, stellar feedback on 
a galactic scale slows down the star formation in the disc 
‘negative’ feedback) ( [Stinson et al.pOO^ hereafter S06) and 


produces outflows of hot gas from the disc to the halo (Dalla 
jVecchia fc Schaye|2008| ). On a cosmological scale, stellar and 
AGN feedback are used to solve the mismatch between the 
halo mass function and the galaxy mass function by quench¬ 
ing star formation in both low-mass and high-mass systems. 
Additionally, stellar feedback pollutes the circumgalactic 
medium (GGlvQ) with metals and produces absorption-line 
systems consistent with those observed around galaxies in 
the spectra of background sources (e.g. Stinson et al.| 20 T 2 ). 
In general, the tuning of feedback in state-of-the-art cosmo¬ 
logical simulations nowadays allows a much better match 
with the observations with respect to the recent past (see 
Vogelsberger et ^|2Q14 Schaye et al.|2015 ). 


Despite the importance of feedback processes, there is 
no general consensus on how they must be implemented in 
numerical simulations. This is due to two reasons. On the 
one hand, these processes occur on scales (both spatial and 
temporal) that are commonly unresolved by the current gen¬ 
eration of simulations. This leads to a plethora of different 
numerical recipes that attempt to approximate the physics 
of feedback on such ‘sub-grid’ scales, often delivering dif¬ 
ferent outcomes despite the similarities of the initial con¬ 
ditions (Scannapieco et al. 2012 ). On the other hand, the 
amount of mechanical and thermal energy deposited into the 
surrounding gas by these processes is very difficult to con¬ 
strain observationally. This implies that the feedback energy 
is treated as a free parameter, which can be tuned ad-hoc 
to reproduce the observations. Even though state-of-the art 
feedback recipes are attempting to address these issues (e.g. 
Keller et al.||2Q14 ), a final solution has still to be found. 


In this paper we investigate how the properties of a 
Milky Way-like system surrounded by its GGM are affected 
by stellar feedback. Our system is produced by the radiative 


^ In this paper the terms CGM and corona have two differ¬ 
ent meanings: the first refers to all the gas that surrounds a 
galaxy disc, the second specifically refers to the hot (~ virial- 
temperature) plasma that settles in equilibrium around the 
galaxy and is built by the hot-mode accretion. 


cooling of a rotating hot gas component (a corona) embed¬ 
ded in a NEW dark matter halo. The gas at the bottom 
of the potential well cools and settles into a disc, eventu¬ 
ally reaching the density required to form stars and pro¬ 
ducing feedback from supernovae and winds. The system 
constituted by the cold, star-forming disc and the hot GGM 
evolves in isolation for 10 Gyr. Thus the properties of the 
final object are affected only by the interplay between these 
two components. This allows us to study the impact of in¬ 
ternal processes (like stellar feedback from the disc and ra¬ 
diative cooling of the corona) on the evolution of a Milky 
Way-like galaxy. The novelties of our work with respect to 
previous studies of simulations of isolated systems are the 
followings: a) our galaxies are surrounded by an extended 
hot corona, which mimics the hot-mode epoch of gas accre¬ 
tion that galaxies like the Milky Way have experienced since 
redshift of ~ 2 (i.e. lookback time of ~ 10 Gyr); and b) we 
present a direct comparison of the gas component with real 
data of the Milky Way by treating the simulated system like 
an observation. 

The paper is structured as follows. In Section [^ we 
present the simulations used in this work. In Section]^ we 
describe the morphology, the mass distribution, the kine¬ 
matics and the gas circulation in these systems. In Section 
^we focus on the gaseous components and mimic observa¬ 
tions of the simulations in order make a direct comparison 
between them and the Milky Way. In Section we discuss 
the results obtained. We present our conclusions in Section 
[ 6 ] 


2 THE SIMULATIONS 


The simulations are run with the N-body + SPH code 
GASOLINE ( Wadsley et al.|2Q04 ). The initial conditions are 
set as in Roskar et al. (2008). We considered a spherical 


NEW dark matter halo with virial radius (r 2 oo) of 200 kpc 
and virial mass of 10^^ M©. Note that these values are per¬ 
fectly compatible with recent estimates of the Galactic dark 
matter halo (e.g. McMillan||2011 ). The halo has an embed¬ 


ded spherical corona of hot gas containing 10 % of the to¬ 
tal mass and following the same density distribution. This 
corona is initially in hydrostatic equilibrium for an adiabatic 
equation of state, and has a spin parameter of A = 0.065 
( Bullock et al.| 20 ()T Maccio et al.|2007 ) with specific angu¬ 
lar momentum oc R, where R is the cylindrical radius. Both 
gas and dark matter components are described by 1 x 10 ® 
particles. Gas particles have initial mass of 1.4 X 10® M© 
and softening length of 50 pc, while dark matter particles 
have softening of ~ 100 pc and masses of either 10® M© or 
3.5 X 10® M©, depending on whether they are inside or out¬ 
side r 2 oo • The initial gas is comprised of a mixture of hydro¬ 
gen and helium, with no metals. 

At the beginning of the simulation radiative cooling is 
switched-on, allowing the central, densest region of the hot 
gas to cool and settle into a disc. Star formation and stellar 
feedback from SNII, SNIa and stellar winds are implemented 
according to the recipes of S06, and the system is followed up 
to t = 10 Gyr. We use a base timestep At = 10 Myr with a 
refinement parameter 77 = 0.175 such that timesteps satisfy 
(5t = At/ 2 ’^ < Tj^eag^ where e is the softening parameter (set 
to 50 pc for baryons and 100 pc for dark matter) and ag is the 
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acceleration a particle experiences. Gas particles also satisfy 
the timestep condition Stgas = r]couranth/[{l-\-a)c-\-f3/Tmax], 
where rjcourant = 0.4, h is the SPH smoothing length, a is 
the shear coefficient, which is set to 1, /3 = 2 is the viscosity 
coefficient and /Xmax is described in Wadsley et al. (2004). 
The SPH smoothing kernel encloses the 32 nearest neigh¬ 
bours. We use an opening angle for gravity of 0 = 0.7. Star 
formation is triggered at a number density n > 0.1cm“^, 
temperature T < 15, 000 K and provided the gas particle 
is part of a converging flow; efficiency of star formation is 
0.05 per dynamical time. Star particles form with an initial 
mass of 1/3 that of the gas particle. Once the mass of a 
gas particle drops below 1/5 of its initial mass, the remain¬ 
ing mass is distributed amongst the nearest neighbors. Star 
particles are assumed to constitute an entire stellar popula¬ 
tion with a Miller-Scalo ( Miller Scalo|[l979 ) initial mass 
function. Feedback from SNII follows a blast-wave recipe 
where thermal energy is injected into the surrounding gas 
and, in addition, gas particles within the blast-wave radius 
have their cooling switched off for a time of the order of 
~ 10^ yr (for details see S06). Stellar feedback also pollutes 
the interstellar medium with metals, whose main impact in 
the simulation is to affect the cooling timescale of the gas. 
We adopt a metallicity-dependent cooling function using the 
prescription of Shen et al. (2010). Metal cooling rates depend 
on density, temperature and metallicity and computed from 
tabulated rates. Metallicity cooling would lead to very low 
temperatures (T ~ 10 K) if allowed to proceed unhindered. 
However below ~ 1000 K the Jeans mass then becomes com¬ 
parable to Mq for reasonable densities, which is much below 
our mass resolution. In order to prevent the artihcial frag¬ 
mentation that would result, we follow Agertz et al. (2009) 
in setting a pressure floor P — . 


As we are interested in distinguishing between gas par¬ 
ticles that have been ejected into the halo by stellar feedback 
and those that remain in the CGM throughout, we do not 
allow metals (or thermal energy) to diffuse from a given gas 
particle to the surrounding ones. With this choice, all par¬ 
ticles have zero metal abundance as long as they are not 
directly affected by SN feedback and stellar winds from the 
disc. In order to test whether this choice severely affects our 
results, we re-ran one of the simulations (F40, see Section 
by implementing metal and thermal diffusion and we found 
no signihcant differences in terms of star formation history 
and star/gas surface density prohle with respect to the non- 
diffusive version. Therefore we preferred to neglect thermal 
and metal diffusion in the simulations. 


The version of GASOLINE used in this work does not 
include the most recent improvements made by |Hopkins| 
(2013) or Hobbs et al. (2013) in the SPH numerical scheme. 
This choice makes it easier to compare our hndings with the 
results obtained in previous studies. Hopkins pointed out 
that, in astrophysical contexts, sub-grid physics is what pri¬ 
marily shapes the outcome of a simulation rather than the 
numerical scheme adopted. Our simulated galaxies do not 
seem to show overcooling in their circumgalactic medium, 
an issue that seems to be caused by the inefficient phase 
mixing of the classical SPH scheme. In addition, the lack 
of signihcant differences between the run with and without 
thermal diffusion indicates that the details of gas mixing are 
not crucial to our setup. 



Figure 1. Comparison between the star formation history of the 
simulated galaxies (thin lines with points) and that of the Milky 
Way determined by|Aumer &: Binney| ( |2009| ) (shaded region) and 
[Fraternali &: Tomassetti] ( |2012| ) (thick grey line). 


2.1 Star formation history and masses 


Using the setup described, we produce four different models: 
F80, F40, FIO and F2.5. Each model uses the same initial 
conditions but differs by the amount of (thermal) energy in¬ 
jected by SNe H into the interstellar medium (ISM), with 
F80 being the most energetic model (80% of the canonical 
10^^ erg per supernova is transferred to the ISM). The main 
physical parameters of the four galaxies after 10 Gyr of evo¬ 
lution are presented in TableNote that, in terms of stellar 
and halo masses, all our simulated galaxies are comparable 
to the Milky Way (e.g. |Sofue et al.||2011 ). 

Fig. [^compares the star formation history (SFH) of all 
the simulations with those derived for the Milky Way by 


Fraternali & Tomassetti 


Aumer & Binney (2009) (double exponential model) and 
( |2012| . The former SFH is based 


on the kinematics and colours of the stars in the solar neigh¬ 
borhood, while the latter is valid for the whole Galaxy and 
is obtained by assuming a linear trend with lookback time. 
Note that the model of Aumer & Binney is normalised to 
the SFR of our s ystems at t=10 Gyr, while |Fraternali 


Tomassetti 


(2012) assume a current SFR of 3M0yr 




In 

all models the SFR is in good agreement with that inferred 
by Aumer & Binney (2009) at all times, and only slightly 


above that predicted by Fraternali & Tomassetti for the last 
~ 7 Gyr. This is an important starting point for the detailed 
comparison between the Galaxy and the simulations. 


2.2 Neutral gas fraction 


Since the simulations did not use radiative transfer, there is 
no direct way to distinguish between neutral and ionised gas. 
In order to estimate the neutral and ionised gas fractions, 
we assume simple collision ionisation equilibrium (GIF), and 
dehne as ‘neutral’ the gas particles with a tem perature be¬ 
low 2 X 10^ K (see Sutherland Dopita||l993 ). The neutral 


hydrogen mass is then derived by assuming a universal hy¬ 
drogen abundance fraction of 0.75. We make no distinction 
between atomic and molecular gas in the simulations. 

A more elaborate approach is the following. Using N- 
body+SPH cosmological simulations post-processed with 


radiative transfer, Rahmati et al. (2013) found that, at 
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Table 1. Properties of the simulated galaxies after lOGyr of evolution. 


Simulation 

F80 

F 40 

FIO 

F2.5 

units 

Mhi+h2 

2.9 X 10^ 

2.5 X 10^ 

2.4 X 10^ 

2.5 X 10^ 

Mq 

Mstars 

5.8 X IQio 

5.9 X IQio 

6.1 X IQio 

6.2 X IQio 

Mq 

SF rate" 

4.2 

4.2 

4.4 

4.5 

Mgyr-l 


® evaluated by considering the amount of stars born in the last 50 Myr. 


each epoch, a relation between particle density and photo¬ 
ionization rate exists. We use this relation (tuned for redshift 
z = 0) to derive the neutral gas fraction for each gas particle 
in the simulations, assuming an UV background radiation 
field based on the model of|Haardt & Madau ( 2001| ) (for 


modify the angular momentum distribution in the corona. 


details see Appendix A of Rahmati et al.||2013 ). The com¬ 
parison between the total maps of neutral gas derived with 
these two approaches revealed almost no differences, thus 
for simplicity we chose to adopt the former simpler method. 


3 RESULTS 
3.1 Morphology 

Fig.|^ shows the final images of the neutral hydrogen 
(atomic+molecular) and of the stellar component, in face- 
on and edge-on projections, for all the simulated galax¬ 
ies analyzed. The maps of the gas component are derived 
by smoothing each particle using a Gaussian kernel with 
FWHM equal to the smoothing length of that particle. 
This gives to these maps a ‘smooth’ appearance that is not 
present in the maps of the stellar component. 

The face-on view reveals that the gaseous discs are ex¬ 
tended roughly as much as their respective stellar compo¬ 
nent, ranging from 12kpc in F2.5 up to 15kpc in F80 at 
the column density level of 10^^ cm“^ (outermost black con¬ 
tour). Sizes do not increase further at 10^® cm“^ (outermost 
grey contour), suggesting that the cold gas discs are trun¬ 
cated and gas at larger radii is too hot to be neutral. This 
truncation is not due to the temperature threshold adopted 
to define the neutral gas phase, as maps derived by using 


a more refined approach (see Section 2.2) show the same 
features at these column densities. Spirals are visible in all 
the gas maps, but the contrast between arm and inter-arm 
regions increases with decreasing feedback. An extreme case 
is F2.5, where the inter-arm regions show extended holes 
in the neutral gas distributions. When seen face-on, stellar 
discs are relatively more similar to each other. Spiral fea¬ 
tures in the gaseous and stellar component of F2.5 are well 
correlated. 

Compared to the size of the HI disc of the M ilky Way 


(30kpc at the column density of 10^® cm , see Kalberla 


Kerp||2QQ9 ), the discs of cold gas in the simulations are 


too small. In general, they appear to fall outside the Hi 
mass-scale relation (e.g. Broeils k, Rhee||1997 ). Since in the 
models the discs of cold gas are produced by the cooling of 
the spinning hot halos, their size will depend on how the 
angular momentum is spatially distributed in these media. 
We have assumed that rotation in the corona is distributed 
as ~ R, which is arbitrary and probably influences the size of 
the gas disc. Additionally, mergers and cold flows may also 


and in fact Wang et al. (2014) successfully reproduced the 
HI mass-scale relation by using cosmological simulations in 
ACDM paradigm. Hence, we do not consider the size of the 
HI discs as a prediction of the simulations. 

The most striking features arise when the systems are 
projected to an edge-on view. Gaseous discs respond to 
the change in feedback by modifying their thickness. This 
happens because high values of Rsn imply a larger injec¬ 
tion of thermal energy inside the disc, which produces more 
turbulence in the ISM and consequently thickens the gas 
layer. A few extra-planar gas features are visible already in 
F40 and FIO at column density lower than 10^^ cm“^; how¬ 
ever, in F80 a thick extra-planar component is visible above 
10 ^^ cm“^ and is even more remarkable at 10^^ cm“^, where 
it extends up to ~ lOkpc from the midplane. This thick 
extra-planar component is remarkably similar to that ob¬ 


served in HI for the nearby disc galaxy NGG891 (Oosterloo 
et al.|2007 ), although its column density is about one order 
of magnitude lower. In Section [T2| we show that not only the 
morphology but also the kinematics are similar. The stellar 
components follow the opposite trend, as stellar discs get 
significantly thicker (i.e., hotter) as feedback decreases. In 
Section [5.II we show that our runs with lower feedback tend 
to form more stars at early times. These stars have a longer 
time to heat vertically, ending up at larger heights. A more 
detailed analysis of this feature will be presented elsewhere. 

3.2 Mass distribution and kinematics 

The mass distribution and the kinematics of the simulated 
galaxies are derived by using an approach similar to the 
‘tilted ring’ method that is often adopted to model the ve¬ 


locity fields of Hi observations (e.g. Begeman 1987). We 
focus on the disc of neutral gas, and divide it into a series 
of concentric rings, each one with a given distance with re¬ 
spect to the centre. This latter is unique and it is given by 
the mass centre of the dark matter distribution. A generic 
ring at distance R from the centre is oriented perpendicular 
to the angular momentum vector of the neutral gas par¬ 
ticles located at that distance from the centre. Therefore, 
the difference of our approach with respect to the classical 
tilted-ring method is that the rings follow the actual three 
dimensional distribution of gas particles rather than the in¬ 
clination and position angle inferred from the velocity field. 
We follow the orientation of each ring in order to properly 
evaluate the rotation velocity, the velocity dispersion and the 
mass surface density for all the components as a function of 
radius. This approach can easily capture symmetric warps 
of the gas component, although non-axisymmetric features 
may not be treated properly. However, none of the models 
show warps in the neutral gas distribution (see Fig.[^. Also, 
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F80 F40 F10 F2.5 



Figure 2. Final images of cold gas and stars of the four simulated galaxies after 10 Gyr of evolution. All maps are on the same scale. 
First row: total maps of cold gas as seen in a face-on projection. Black contours are at column densities of 10^®, 10^°, 10^^, 10^^ cm“^; 
grey outermost contours are at 10^® and cm“^. Second row: maps of the stellar component as seen in a face-on projection. Third 

and fourth rows: as the first and second rows, but the galaxies are seen in a edge-on projection. 


since this approach relies on the neutral gas alone to infer 
the inclination of the various rings, it can fail if the disc of 
stars and gas are significantly misaligned. This is also not 
the case for these models. 

Fig.[^ shows the gas and stellar surface densities (top 
row) and the decomposed rotation curves of cold gas (bot¬ 
tom row), and a comparison with a mass-decomposition of 
our Galaxy. To derive this latter, we assumed that the Galac¬ 
tic gravitational potential is produced by four components: 
a stellar bulge, an double-exponential stellar disc, a NFW 
dark matter halo and a disc of cold gas. Despite the strong 


evidence for a ‘peanut-shaped’ boxy-bulge in the Milky Way 


(e.g. Dwek et al.|1995 Nataf et al.|2QlQ McWilliam & Zoc- 


cali|2Q10' Ness et al.|2012 ), a de Vaucouleurs bulge provides 

a very good fit to the Galaxy’s rotation curve (Sofue et al. 


2009). As modelling the stellar distribution and kinematics 
of the inner Galaxy is beyond the purpose of this work, for 
simplicity we decided to adopt the classical de Vaucouleurs 
profile to describe the Galactic bulge. We fitted the param¬ 
eters of these components to the Milky Way’s rotation ve¬ 


locity data of Sofue et al. (2009), which we rescaled to more 
up-to-date values of the Galactic constants i ’0 = 239kms“^ 
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Figure 3. First row: surface density profiles of the various components for the simulated galaxies (first four columns), compared with 
that of the Milky Way (last column). Second row: Hi rotation curve and contribution of the various components to the circular velocity 
for the same galaxies. The Milky Way’s rotation curve is t aken fro mlSofu e et al. | ( |2009| > but it has been rescaled for the following values 
of the Galactic constants: vq =239kms“^, Rq =8.3kpc ( McMillan||2Qll >. The cold gas (H 1 +H 2 , corrected for He) surface density of 
the Milky Way is taken from |Binney &: Merrifield] ( |1998| >, while the stellar surface density and the circular velocities are derived from 
our mass decomposition (see text for details). 


and Rq — 8.3kpc (McMillan 2011). As the rotation curve 


in the outer disc is very uncertain, we assumed that the ro¬ 
tation velocity flattens at i? > Rq and discarded the data 
points in that region (light-grey points in Fig.|^. We used 
the cold gas (HI+H2, plus a correction for the He) surface 
density of |Binney fc Merrifield (1998) and we kept it fixed 
in our fit. In addition, we constrained the total baryonic sur¬ 
face density at R — Rq to be in the range 54.2 zb 4.9 M© pc“^ 
( Read|20T4 ). The details of this mass decomposition will be 
presented in Marasco & Fraternali (in prep). 

The cold gas (HI+H2) surface density profiles of the 
simulated systems are all flat at the level of ~ IOMqpc”^ 
(or ~ 10^^ cm“^), and increase by one order of magnitude 
close to the centre where star-forming gas piles up. All pro¬ 
files are truncated, as at about 10 kpc the surface density 
decreases by two orders of magnitude in a small (2-3 kpc) 
region. This is different from what has been found in the 
Milky Way, where the cold gas fades gently as the radius 
increases. All simulated stellar profiles show the presence of 
a mass concentration in the centre due to a bulge. We notice 
that the central stellar surface density increases slightly with 
feedback, from 6 x 10^ M© pc“^ in F2.5 to 8 x 10^ M© pc“^ in 
F80. Outside the region of the bulge, stars distribute in a disc 
that seems to show, in ah simulated systems, a double expo¬ 
nential prohle, with the break between the two slopes occur¬ 
ring approximately where gaseous discs are truncated. This 
is not surprising: stars located beyond the gaseous discs must 
have been transported to that location by radial migration 


via transient spiral corotation capture (Sehwood & Binney 
2002 Roskar et al. poo's ). A break in the slope of stellar pro¬ 
hle has been observed in several galaxies (e.g. Kregel et al. 
2002) and it is possibly present in the Milky Way as well 
^ Sale et al.|2010 Minniti et al.|20lT ). Note that the Galactic 
Hi disc also seems to be truncated at R ~ 15kpc (Kalberla 


& Dedes 2008), although this truncation is shallower than 


that shown by our simulated galaxies. Truncated HI discs 
are sometimes observed in external galaxies, but they are 
thought to be caused by photo-ionization from the extra- 


galactic UV background (e.g. Maloney|1993 ) rather then by 
an effective drop in cold gas density. 

The rotation curves, derived by measuring the actual 
rotation velocity of cold gas particles in the various rings 
(black lines), show several interesting features. Ah rotation 
prohles are steeply rising up to 350 — 400 km and then de¬ 

cline and hatten to about 250kms“^ beyond R^3kpc. The 
Milky Way rotation curve has a similar shape, but the inner¬ 
most peak reaches a lower value (slightly above 250kms“^) 
and it battens to a slightly lower rotational speed, whose 
precise value depends on the choice of the local standard 
of rest rotational velocity (here assumed to be 239kms“^, 
see McMillan 2011). Using the particle distribution in the 


simulation it is possible to compute the gravitational po¬ 
tential and therefore to derive the circular rotation for ah 
the mass components separately. In ah models, the dynam¬ 
ics of the system is driven by the stellar component out to 
the radius where the gaseous disc is truncated. Beyond that 
point, the dark matter becomes the dominant component, 
while gaseous discs are not dynamically dominant at any ra¬ 
dius as is expected in this type of galaxies. As shown in the 
lower panels of Fig.[^ the neutral gas is in perfect circular 
rotation except in the outermost regions of F80. However, as 
we show is Section [4.2| part of the cold gas in these regions 
is extra-planar and rotates at a slower speed than the gas in 
the midplane. 

The shape of the dark matter prohle deserves mention. 
While the outermost part is similar in ah models, the in¬ 
nermost prohle rises signihcantly more steeply when higher 
feedback energy is coupled to the ISM. This corresponds to 
a higher mass concentration, baryonic and non-baryonic, in 
the centre of galaxies with high feedback, and to a higher 
peak in rotational velocity. This result is surprising as it 
is opposite to what we expected. Higher feedback has been 
historically invoked to wash out the mass built-up in the 


centre of systems (the so-called core-cusp problem, e.g. de 


Blok 2010): as baryonic matter is removed from the sys¬ 
tem centre by stellar feedback, dark matter follows it (e.g. 
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Figure 4. Top panel: velocity dispersion profiles for both stars 
and cold gas in all the models. Bottom panel: scale-height of the 
cold gas as a function of radius in the simulations. 


Navarro et al.|p^^ [Read &: GilmQre||20Q5| |Ogiya &: Mori| 
201 1| ). Our simulations show the opposite trend, as the cen¬ 
tral stellar surface densities increases with and conse¬ 
quently compressing the dark matter more in the systems 
with larger feedback output. This indicates that the feed¬ 
back scheme adopted is not very efficient in removing low- 
angular momentum baryons from the system centres: in fact, 
it promotes their growth. We discuss the impact of stellar 
feedback on our simulations in more detail in Section [5T] 
The top panel of Fig.|^ shows the stellar and neutral 
gas velocity dispersion profile, cr* and (Jgas, for the different 
models, computed as 


= \j\ [c7l{R) + am) + ^l{R)) 

= \j\ {al{R)+al{R) + ol{R)) 


KbT 

mu 


( 1 ) 

( 2 ) 


where an, Gz and Gcf are particle-to-particle velocity dis¬ 
persions in the various directions (radial, vertical and az¬ 
imuthal), averaged over all neutral gas particles in each ring, 
T is the mass-averaged temperature of these particles, and 
ttih is the hydrogen mass. 

In all models, cTgas is flat and its value depends on the 
feedback used, with F2.5 being the least turbulent system. 
F80 has very turbulent motions so that cTgas increases up to 
25 — 30kms“^, ~ 2 — 3x larger than that measured for the 
Hi in our Galaxy and in nearby disc galaxies with similar 
star formation rates (e.g. Boomsma et al.|2008 ). The veloc¬ 
ity dispersion of the stellar component is larger than that of 
the gas and it increases in the inner regions where the bulge 
dominates the total surface density. As the edge-on maps 
suggest (bottom panels of Fig.[^, the thickening of the stel¬ 
lar discs in the models with lower feedback corresponds to 
a general increase in cr*. At 15kpc, there is a factor of 
~ 2.5 of difference in a* between F2.5 and F80. 

The bottom panel of Fig.|^ shows the scale-height of 


cold gas as a function of R for the various models. All 
simulated galaxies show a significant flaring, which is more 
prominent in the models with higher feedback. The HI scale- 
height in the inner disc of the Milky Way is between 100 and 
200 pc ( Dickey Lockm^|1990 ), so the cold gaseous disc 
of F80 is thicker than that of our Galaxy. However, we point 
out that the gas velocity dispersion, and therefore the scale- 
height, are resolution dependent. These two quantities are 
expected to correlate for discs in hydrostatic equilibrium, 
and in SPH simulations gravity - and therefore hydrostatic 
equilibrium - cannot be described accurately on scales of the 
order of the gravitational softening e. In all the simulations 
we adopted e = 50 pc for the stars and gas, implying that 
gravity deviates from Newtonian on scales comparable to 
the scale-heights of HI discs in galaxies. Thus we stress that, 
even though the trends shown by the scale-height and veloc¬ 
ity dispersion with feedback in the models are certainly gen¬ 
uine, the actual values of these quantities are probably over¬ 
predicted as the gravitational restoring force to the plane is 
underpredicted within a few e from the midplane. 


3.3 Time-evolution of stars and neutral gas 

We evolve the models for lOGyr. Fig.|^ shows the evolu¬ 
tion of the systems. Thinner lines represent earlier stages 
of the evolution (i.e., larger lookback times). In all systems 
both gaseous and stellar discs form inside-out. Also, at all 
times, stellar discs are systematically more extended than 
cold gaseous discs, with the exception of F80 and F40 at 
t = 1 Gyr. Stars are born within the cold gas disc then ra¬ 
dial migration drives stars outwards, producing the system¬ 
atic mismatch in sizes. We already noticed in Section |3.2| 
that stellar discs show a double exponential profile, with 
the break in the slope occurring at the point where gaseous 
discs are truncated. Fig.shows that this feature occurs at 
any time. Stellar migration seems to continuously transfer 
stars from the place where they are born, the gaseous discs, 
to outer radii, thus producing the observed double exponen¬ 
tial distribution in the stellar profile (see |Roskar et al.|2008 ). 
Along with stellar discs, gaseous discs also grow inside-out, 
extending to larger radii at later times. This happens in the 
simulations because the gas cooling time increases outwards 
where the density is lower, thus cold gas at large radii be¬ 
comes available only during the later stage of the evolution 
of the system ( Roskar et al.|2Q0'8 ). 


4 DIRECT COMPARISON WITH THE GAS OF 
THE MILKY WAY 

In this Section, we focus on the gas component. Our ap¬ 
proach is to simulate an ‘observation’ by considering an ob¬ 
server placed inside the disc of the simulated galaxies at 
R = 8.5 kpc from the centre, and studying how - from this 
specific point of view - gas particles distribute as a func¬ 
tion of the angular position in the simulated sky and of the 
line-of-sight velocity in the local standard of rest (LSR). 
This ‘pseudo’ LSR is determined by averaging the motion of 
about 100 star particles around the selected point of view. 
We stress that this approach allows a direct comparison be¬ 
tween the gas component of the simulated systems and that 
































Marasco et al. 



Figure 5. Time-evolution of stellar surface density and cold gas surface density in the simulations. Thicker lines represent later stages 
of evolution, as indicated in the leftmost panel. 


of our Galaxy, for which only 3D position-velocity informa¬ 
tion are available. 

We compare this simulated observation with two dif¬ 
ferent datasets. The neutral gas (log(T) < 4.3) is com¬ 
pared with the all-sky Hi emission data of the LAB Sur¬ 
vey ( [Kalberla et al.|[20Q5| , whereas gas at higher tempera¬ 
tures (4.3 < log(T) < 5.7) is compared with the absorption 
line measurements of ionized species of Lehner et al. (2012), 
Sembach et al.| ( [20Q3| and [Savage et al. (2003). 


4.1 The Hi disc 

The LAB Survey is an all-sky datacube of Galactic (IulsrI < 
400kms“^) Hi emission. In order to compare the simula¬ 
tions with the LAB data, we produce a simulated datacube 
(a ‘modelcube’) by using the following procedure: 

• we evaluate the pseudo-LSR as described above; 

• in this reference frame, we evaluate the longitude /, the 
latitude b and the line-of-sight velocity ulsr of each parti¬ 
cle. The latitude 6 = 0° is chosen to be aligned to the disc 
midplane, so that this reference frame effectively mimics the 
Galactic coordinate system. 

• We place each particle in a 3D grid (/,6, u) with pixel 
size equal to that of the LAB data; 

• we smooth the neutral gas content of each particle 
in velocity according to its thermal velocity dispersion 
{kBT//j./mp)^. 

• we smooth the neutral gas content of each particle ac¬ 
cording to the desired spatial resolution (see below); 

• we convert each pixel of this modelcube to a brightness 
temperature by assuming an optically thin regime. 

The spatial smoothing process deserves further discussion. 
The ‘natural’ resolution of the simulations is the SPH ker¬ 
nel size, that is the distance between a given particle and its 
32^^^ neighbour. If we smooth each particle separately to a 
spatial resolution equal its SPH kernel size, we assure that 
the whole simulated box is sampled without discontinuities. 
However, the smoothing lengths vary significantly from the 
galaxy centers, where particles crowd and the resolution is 
high, to the outskirt of the systems, where particles are rarer 
and the resolution is significantly lower. This effect is dra¬ 
matically amplified by the inner projection used: we found 
that, on the midplane, the angular resolution varies from less 
than 2° at / = 0° (i.e. in the direction of the galactic centre) 
up to ~ 15° at / = 180°. In addition, the angular resolution 


gets significantly worse at higher latitudes. As we are com¬ 
paring our modelcubes with data that have a fixed angular 
resolution, we decided not to use this smoothing method and 
instead to smooth each particle to a fixed angular resolution 
by using a 2D Gaussian kernel with a FWHM of 8° x 8°. 

Fig.|6] shows longitude-velocity {I — v) slices taken at 
three different latitudes (±15° and 0°) for our modelcubes, 
and compare them with the LAB data smoothed at the same 
spatial resolution (8°). The HI line profiles of all cubes 
have been smoothed in velocity by using a Gaussian ker¬ 
nel with a FWHM of 10kms“^ in order to emphasize the 
low-level emission in the LAB data. The midplane (6 = 0°) 
panels show that the velocity spread of the HI emission de¬ 
creases with decreasing feedback, which can be explained by 
the drop in cTgas discussed in Section [3^ The main differ¬ 
ences between the simulated galaxies and the Milky Way 
at this latitude are that HI emission in the former reaches 
larger velocities around I = 0° and lower velocities around 
/ = 90°, 270°. This happens respectively because a) in the 
simulations the rotation curves peak at higher velocities 
than those reached in our Galaxy; b) the gaseous discs are 
smaller than that of the Milky Way. These considerations 
apart, from the slice at 6 = 0° alone it is not possible to de¬ 
cide which model better resembles the Galactic HI emission. 

Moving to different latitudes (top and bottom rows of 
Fig.§ it becomes clear that, as feedback decreases, the l—v 
plots lose more of the sinusoidal shape that is visible in the 
data. This shape is produced by the rotation of the HI layers 
located immediately above the midplane and results from 
the vertical extent and the kinematics of this gas. In the 
Milky Way, most of the low-level emission visible at these 
latitudes is due to the extra-planar H1. This medium, which 


is thought to be produced by the Galactic fountain (Marasco 
et al.||2Q12 ), rotates slower than the gas in the disc and ex¬ 


tends up to a few kpc above the midplane ( Marasco Fra- 
ternalipOll ). Both these features significantly enhance the 


sinusoidal shape visible in the LAB data. As a comparison, 
on each panel we overlaid a thick red contour representing 
the HI emission (at the level of 0.01 K) predicted for a dif¬ 
ferentially rotating disc extended up to 25 kpc and with a 
scale-height of 150 pc, approximately the value estimated 


for the thin HI disc of our Galaxy out to R = Rq (Dickey & 
Lockman|1990 ). For simplicity, we do not attempt to model 


the warp in the outer HI Galactic disc (e.g. Levine et al. 


2006), which at any rate would affect only a small portion of 
our modelcube (40°</<160°, —200 < ulsr <—50 km s“^). 
Glearly, most of the low-level emission of the LAB data is lo- 
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Figure 6. Longitude-velocity slices taken at three different latitudes (as indicated on top of the leftmost panels) for cold gas in the 
simulated galaxies (columns 1 to 4) and the LAB data of the Milky Way (5th column). Black contours are at brightness temperatures 
ranging from 0.01 K to 40.96 K in multiples of 4. In the LAB data an additional contour at —O.OIK is shown in grey. The red contour 
represents the HI emission at the level of O.OIK from a model of Galactic disc with a scale-height of 150 pc. All panels have been 
smoothed to 8° of angular resolution and 10kms“^ of velocity resolution. 


Gated outside this contour. Neither F2.5 nor FIO have extra- 
planar gas (Section |3.1| and their l — v emission is confined 
within the contour of the thin disc model. F40 also does not 
have extra-planar gas, but the average thickness of its disc is 
larger than 150 pc (bottom panel of Fig.[^ thus some emis¬ 
sion leaks outside the red contour. Finally, F80 has a layer 
of extra-planar gas. Even though at this stage we do not 
know what the kinematics of this layer is, the fact that l—v 
emission of F80 is the one that resembles the LAB data best 
would suggest that the extra-planar material rotates slower 
than the gas in the thin disc. We show that this is the case 
in the next Section. 


4.2 The extra-planar HI 

In the last 10-15 years, deep HI observations reaching col¬ 
umn densities of ~ 10^^ cm“^ have proved that star-forming 
disc galaxies often show a vertically extended extra-planar 


Hi component (or Thick HI discs’, e.g. Fraternali et al. 2002 


Barbieri et al.|20'05 Gentile et aL|2013 ). Our Galaxy is not 

an exception ( [Marasco fc Fraternali|2Qll ). A detailed mod¬ 
elling of this component reveals that its rotational velocity 
decreases with increasing distance from the midplane. These 
kinematics can be explained by a model of the galactic foun¬ 


tain (Bregman 1980) interacting with an ambient medium 


at a lower specific angular momentum than that of the disc 
( Fraternali Binney||2QQ8 ). 

F80 is the only simulated system that shows a promi¬ 
nent thick HI disc. We directly evaluate the rotation curves 
above the disc: Fig.[^ compares the rotation curve of cold 
gas derived in the midplane of F80 with those obtained at 
larger heights. Rings with less than 10 cold gas particles have 
been neglected. The figure clearly shows that the larger is 


the distance from the midplane, the slower is the rotation. 
At 2kpc above the midplane, the difference in rotational 
speed with respect to the midplane is about 30kms“^, cor¬ 
responding to a lag of ~ 15kms“^ kpc“^, in remarkable 
agreement with that found by Marasco & Fraternali (2011) 
for the HI halo of the Milky Way and by |Oosterloo et aL| 
(2007) for NGG 891. We discuss the origin of this rotational 
lag in Section |5.2| Fig.[^ also shows that this component 
is not present in the innermost region of the disc, point¬ 
ing to an HI halo with a toroidal shape. This is expected in 
the galactic fountain mechanism, as the gravitational restor¬ 
ing force of the disc decreases outwards making it easier for 


HI gas to reach larger distances from the plane (Fraternali 


& Binney 2006). We found that all the cold gas particles 


that form the extra-planar gas layer have non-zero metallic- 
ity, indicating that this material originated inside the disc. 
Specifically, the metal abundance of this component is the 
same as the underlying disc at all radii, which implies that 
neutral gas ejected from the disc by stellar feedback does 
not significantly change its galactocentric distance during 
its trajectory. 


In our systems, cold gas is absent at large distances from 
the discs (see Fig.|^. On the one hand, this is in contrast 


with the results of the GOS-Halo Survey (Tumlinson et al. 


2013 ), which revealed that cold, photo-ionized absorption 
systems are widespread around galaxies ( Werk et al.||20I3' 


20I4| ). On the other hand, the dense HI clumps that were 


ubiquitous in the GGM of galaxies in previous SPH simula¬ 
tions (e.g. Kaufmann et al. 2006 2009) are not present in 
our simulations. This in agreement with the results of [Pisano] 


et al. (2011), who show the lack of massive HI floating clouds 


at large distances from galaxies in Local-group-like systems. 


























































10 Marasco et al. 


F80 - kinematics of the extraplanar gas 



Figure 7. Rotation velocity of neutral gas in model F80 at dif¬ 
ferent heights above the midplane. 


4.3 The warm+hot gas 


By analyzing the spectra of background sources, a num¬ 
ber of studies have revealed the presence of absorption sys¬ 
tems of ionized species (Silll, SilV, CiV, OVI) at |i;l srI < 
400kms“^ (e.g. Wakker et al.|20d3 Lehner et al. 2012). The 
global kinematics of these absorbers is not consistent with 
being part of the Galaxy disc, which makes of these sys¬ 
tems an ionized counterpart of the classical HI High-Velocity 
Clouds. All these ionized species are expected to exist at a 
temperature well below the Milky Way’s virial temperature. 
Hence, a natural explanation for these absorbers is that they 
are produced by the cooling of the CGM of the Milky Way 
and represent a strong evidence for accretion of ionized gas 
onto the Galaxy disc ( Lehner V Howk|2011 Fraternali et al. 


2013). 


We now compare the position-velocity distribution and 
the column densities of these systems with those predicted 
for the warm-hot (4.3 < log(T) < 5.7) gas in the simulation 
with the highest feedback (F80) which already shows a layer 
of extra-planar cold gas in good agreement with the obser¬ 
vations. In addition, we use the full 6D information available 
in the simulation to determine the location of the absorbing 
material. As we will remark during our analysis, our results 
remain valid also for the runs with lower feedback. 


4 . 3.1 Kinematics 


The analysis that we perform here is similar to that used 
by Marasco et al. (2013) to compare their galactic fountain 
model with the absorption data. In the following, we summa¬ 
rize the procedure adopted. First, we use the same method 
described in Section [T^ to produce modelcubes of the warm 


(4.3 < log(T) < 5.3) and hot (5.3 < log(T) < 5.7) gas of F80. 
The main difference is that here the spatial smoothing uses 
the SPH kernel size of each particle instead of a fixed an¬ 
gular resolution: as our data consist of a set of absorption 
line measurements that are randomly distributed in the sky, 
we prefer to sample the warm-hot gas in the simulated box 
smoothly and extend its ‘filling-factor’ as much as possible. 
Each modelcube is produced twice: once using only the pol¬ 
luted gas (i.e., with metal abundance larger than 0), and the 


second time using both polluted and unpolluted gas. We re¬ 
mind the reader that, as metal diffusion is not used and the 
simulation starts with no metals, particles can be polluted 
only by SN explosions in the disc or stellar winds. There¬ 
fore, the first modelcube is representative for gas that has 
been affected by feedback and participates in the fountain- 
cycle. A qualitative comparison between models and data is 
carried out by overlapping the position-velocity of the ab¬ 
sorbers on top of the I — v diagram of the modelcubes. A 
more quantitative analysis uses an iterative KS-test to de¬ 
rive the fraction of detections that is consistent with our 
modelcubes (for the details see Marasco et al. p0l3l l. 

Fig.|^shows two representative longitude-velocity {l—v) 
slices for the warm (panels on the left) and the hot (panels 
on the right) gas of F80. Slices derived for the polluted gas 
alone and for the polluted+pristine material are adjacent to 
each other. The three contours shown enclose 68.3%, 95.4% 
and 99.7% of the total flux of the modelcubes, by analogy 
with a Gaussian distribution. The location of the observed 
absorption systems in the position-velocity diagram for these 
latitudes is plotted on top of the slices. As in [Marasco et al.| 
(2013), the modelcube of warm gas is compared to the ab¬ 
sorption data of Lehner et al. ( 2012| ) (Si III, Si IV, GlI, Gill, 
GIV) while the hot gas is compared to the joined datasets 
of Sembach et al. ( 2003| and Savage et al. ( 20Q3| (OVi). 
Warm material at Ir’Lsnl <90kms“^ has been neglected as 
it was excluded in the observations of Lehner et al.j ( |2012| . 

The latitude bins chosen emphasize the differences be¬ 
tween the whole gas and the polluted material alone. Sev¬ 
eral detections occur at velocities that are well beyond those 
predicted by the warm-hot polluted material of the simula¬ 
tion. Instead, the inclusion of pristine gas from the GGM, 
which rotates at a lower speed than the disc, significantly 
increases the relative line-of-sight velocities reached by the 
warm-hot material and produces a better overlap with the 
observations. Using our iterative KS test we find that only 
6% of the detections of Lehner et al. are compatible with 
the warm polluted material alone, whereas this percentage 
increases to 45% when including the unpolluted gas. This 
suggests that the gas cycle produced by stellar feedback can 
not be solely responsible for the observed warm absorption 
features, as many absorbers have kinematics that are more 
consistent with that of the GGM. This discrepancy is signif¬ 
icantly mitigated when we consider the OVI detections of 
Sembach and Savage, as the percentage of reproduced fea¬ 
tures passes from 31% for the hot polluted material alone to 
51% when all the gas is considered, indicating that a good 
fraction of OVI around galaxies might come from stellar 
feedback. We notice that these results do not change when 
metal and thermal diffusions are included in the simulations: 
we derived the same l — v slices of Fig. [fusing the version 
of F40 where diffusion is included, and we found similar re¬ 
sults. Still, the models can account at most for about half of 
the detections. A possibility is that some of the high-velocity 
O VI absorption features observed by Sembach et al. (2003) 
could originate in the GGM of nearby galaxies and are there¬ 
fore not related to the Milky Way. Following [Marasco et ah] 
(2013), we removed from our catalogue 41 O VI features that 
can be associated to external galaxies, and found that the 
fraction of the remaining O VI detections that is reproduced 
by our model increases to 70%. We therefore conclude that 
the kinematics of the hot gas predicted by our runs is com- 
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warm gas (polluted only) 


warm gas (all) 


hot gas (polluted only) 


hot gas (all) 



Figure 8. Longitude-velocity slices in two different representative latitude bins (indicated on the top-right corner of each panel) for the 
warm (4.3 < log{T) < 5.3, left-hand side panels) and hot (5.3 < log{T) < 5.7, right-hand side panels) gas of our run F80 (filled region). 
Circles overlapped to the warm gas represents the detections from |Lehner et aT] ( |2012| >, those overlapped to the hot gas are taken from 
[Sembach et al.| ( |2003| > and [Savage et ar] ( |2003| >. Contour levels from the innermost to the outermost are at 68.3%, 95.4% and 99.7% of 
the total flux. The HI emission of the classical HVCs is reported in the various panels as the thick contours, labelled with the name of 
the respective complexes. 
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Figure 9. Gas density profile (solid line, left-hand axis) and tem¬ 
perature profile (long-dashed line, right-hand axis) in our model 
F80. The top panel shows the profiles derived radially, i.e. inter¬ 
secting the galaxy disc, while the bottom panel shows those de¬ 
rived vertically, i.e. along the rotation axis. The density profiles 
are decomposed in four components at different temperatures, as 
indicated in the top panel. 


patible with that shown by the majority of the OVI ab¬ 
sorbers observed around the Galaxy. 


4-3.2 Location 

Under the hypothesis that the gas responsible for the ob¬ 
served absorption features is best represented by a mixture 
of warm-hot polluted and unpolluted gas, we use the simu¬ 
lation to determine what is the typical distance of the inter¬ 
vening absorption material. 

Fig. shows both the radial (i.e., parallel to the disc) 
and the vertical (perpendicular to the disc) gas volume den¬ 


sity profiles in model F80. The former has been derived by 
considering all gas particles in a cone centred on the system 
centre with an aperture of 5° above and below the midplane, 
which allows us to intersect the region of the disc and ex¬ 
clude any contamination from gas above it, while for the 
latter we used a cone of 30° around the rotation axis of the 
system. We used different lines and colours to represent the 
various phases of the gas. Cold gas is present only inside the 
galaxy disc, as it appears only in the inner ~ 10 kpc of the 
radial profile. The vertical profile does not show the pres¬ 
ence of cold material even in the innermost regions, where 
cold extra-planar gas is located. This happens because, as 
discussed in Section [4.2| this layer has a toroidal shape and 
therefore it is missed by the cone of view that we adopt. 
Warm and hot gas are instead located both in the disc and 
in the CGM. 


As the kinematics of the warm gas coming from the disc 
is not in agreement with the observations (Section [4.3.1 ) we 
must focus on the warm material in the CGM, which appears 
to be located at distances between 150 kpc and 500 kpc from 
the centre, beyond the hot gas. As we are adopting a simple 
temperature threshold to discriminate between warm and 
hot material, this can be interpreted in terms of a temper¬ 
ature gradient in the CGM. The temperature of the CGM 
in the inner region of the system is larger than that in the 
outskirts (long-dashed line in the bottom panel of Fig. © , in¬ 
ducing a spatial separation between warm and hot material 
and creating a deficit of warm gas between 20 and 200 kpc. 
However, gas at this temperature is commonly seen in ab¬ 


sorption within 150 kpc from galaxies (e.g. Tumlinson et al. 


2013). Furthermore, in the Milky Way, there is evidence that 


most of the warm absorbers are located in the lower halo. 


within a few kpc from the Galaxy disc (Lehner & Howk 


2011 Wakker et al.||2012 Marasco et al.||2Q13). 

The hot gas in F80 extends radially from the system 
centre up to more than ~ 300 kpc from the disc, although 
it is absent in the inner 40 kpc in the direction perpendic¬ 
ular to the midplane. This is in agreement with the OVI 
impact-parameters derived in external galaxies: both jStockej 
et al. (2006) and Wakker & Savage (2009) found that low- 
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redshift intergalactic OVI originates within 500 kpc from 
bright galaxies. In order to verify whether the absorption 
systems are widespread in the CGM or are confined within 
a shell at a given distance from the system centre, we pro¬ 
ceed as follows: we consider the angular position and the 


^LSR of the O VI absorbers found by Sembach et al. (2003) 


and Savage et al. (2003), and we determine at what dis¬ 


tance from the disc the hot gas particles of F80 with the 
same position and r’LSR are located. We found that there is 
no preferential location for these particles, ~ 70% of them 
lie between 20 and 200 kpc from the centre. 

Our findings suggest that the location of the hot gas pre¬ 
dicted by model F80 is in agreement with the absorption-line 
observations, while that of the material at a lower temper¬ 
ature is not. We stress that our results are valid also for 
the models with lower feedback, as their density and tem¬ 
perature profiles are very similar to those of F80. The only 
exceptions are found in the inner region of the radial pro- 
hle where, as the feedback decreases, the hot material inside 
the disc gets increasingly replaced by warm gas in response 
to the lower amount of energy injected into the ISM. One 
may wonder whether our findings are affected by the lack 
of metal and thermal diffusion, as the mixing between the 
different gas phases can alter the temperature profile of the 
CGM. To address this question, we compared the gas tem¬ 
perature and density distribution in F40 with those derived 
by re-simulating the same system with thermal and metal 
diffusion included. We found no differences, except that the 
low-density hot gas located between 20 and 80 kpc from the 
center disappears, likely because it gets thermalized by the 
(hotter and denser) surrounding coronal gas. It is possible 
that warm or cold material arises in the inner hundreds of 
kpc from the discs only in a full cosmological context as a 
consequence of late-time filamentary cold-mode accretion, 
as suggested by high-resolution cosmological simulations of 


Milky Way-like systems (Keres & Hernquist 2009 Joung 


et al.|20l^ Stinson et al.|2012 ). Insofar as our findings point 


towards this scenario, we remind the reader that the evolu¬ 
tion of the CGM in the systems depends in principle on the 
initial conditions of the simulation. Adopting an isothermal 
corona (e.g. Binney et al.|[2009 ), a lower baryonic fraction 
(motivated by X-ray observations, see Anderson & Bregman 
[Mot , a different distribution of the initial angular momen¬ 
tum, the presence of substructure, or a different feedback 
recipe are all choices that might have lead to different re¬ 
sults. 


In the models, the OVI mass movi associated to each 
gas particle of total mass m can be written as 


movi = 


/ movi \ 

f TT^heavy ^ 

f mo \ 

V J 

\ m J 

\mheavy J 


(3) 


where mo is the oxygen mass associated to the particle and 
TTiheavy IS thc mass of thc clcments heavier than helium. 
The three bracketed terms in eq. it can be determined sep¬ 
arately. We assume collisional ionization equilibrium (GIF) 
to evaluate (mo^^/mo) as a function of the temperature 
of the particle, by using the oxygen ionization-balance ta¬ 
ble of Sutherland & Dopita (1993). Note that the assump¬ 


tion of CIE should be solid for OVI systems in the halo of 


galaxies, as discussed by Sembach et al. (2003). The quan¬ 
tity {mheavy!'^) IS thc mctal abundance of the particle, 
which is known in all our simulations. Finally, we assume 
solar abundance ratios for our gas and set (mp/ mhea yv) = 
(mo/?TiH) 0 X (mheavy/mii)Q^ = X (0.0191)“^ (Lod- 

|ders||2010| ). 

We use a procedure similar to that described in Section 
4.1 to build a modelcube where each voxel /(/, 5, v) contains 


information on the OVI column density per unit velocity 
at the position (/,5) of the simulated sky. The OVI column 
density at position (/, b) will be simply given by 


No, 


nv 

,{l,b)= / 

J V^ 


I (1, 5, v) dv 


(4) 


where vi and V 2 are two generic line-of-sight velocities. Fol¬ 
lowing Marasco^ta^ ([201^ , we use eq. to compute iVovi 
in all line-of-sights {l,b) where OVI detections have been 
found in the combined datasets of [Sembach et ah (2003) 
and Savage et al. ( 2003| ). The integral in eq. (|4J is calculated 
between v — ^bw and v-\- ^bw, where v is the mean veloc¬ 
ity and bw is the line-width of the observed O VI absorption 
feature. By considering all the OVI column densities com¬ 
puted in the various line-of-sights, we find for model F80 an 
average value of (logA^Ovi) = 13.7 ± 1.03, compatible with 
the value derived from the observations, 14.16 ± 0.34. This 
result changes little when metal and thermal diffusions are 
included, as the average O VI column density in the version 
of F40 that includes diffusion is 14.48 ± 1.45. Gonsidering all 
the uncertainties, we conclude that the systems show a layer 
of hot gas whose kinematics, location and column densities 
are compatible with what is inferred for the Milky Way and 
external galaxies from O VI measurements. 


4-3.3 Column density 

Since the hot phase of the gas in the simulation is the one 
that best matches the position and the kinematics of the 
observed OVI absorption systems, it is interesting to test 
whether the amount of hot gas in the simulations is con¬ 
sistent with the observed OVI column densities. The com¬ 
parison between the simulations and the observations is not 
straightforward, as the former are not designed to keep track 
of the various heavy elements nor to determine the gas ion¬ 
ization balance. In the following, we attempt to determine 
the O VI column density in a generic line of sight A^Ovi(^? 
by using the information that is available for each particle, 
namely its total mass, temperature and metal abundance. 


5 DISCUSSION 

5.1 The impact of stellar feedback 

One of the most puzzling features of the simulated galaxies 
is that, as shown by Fig.^ all star formation histories are 
very similar despite the fact that Fsn changes significantly 
from simulation to simulation. This result dramatically dif¬ 
fers from that found by S06, i.e. that the injection of en¬ 
ergy into the ISM slows down the star formation processes 
and therefore the same amount of fuel gets smoothly con¬ 
sumed on longer timescales. The main difference between 
our work and that of S06 is that the systems have coro- 
nae and hot-mode accretion happens all the time. Our find¬ 
ings indicate that, at least in the range of feedback energies 
studied in this work, gas accretes onto the discs - and star 
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Figure 10. Total stellar mass assembled in our runs, normalized 
to the stellar mass assembled by F80, as a function of time. The 
x-axis is in logarithmic scale to emphasize the early stage of the 
evolution. 


formation proceeds - in a similar fashion for all the simu¬ 
lated systems. But is stellar feedback affecting the hot-mode 
accretion at all? To answer this question, we ran another 
simulation (FO) adopting the same initial conditions as the 
other simulations, but with stellar feedback switched off al¬ 
together. Fig.p^ shows the total stellar mass assembled in 
all our runs (including FO), normalized to the stellar mass 
assembled by F80, as a function of time. Clearly, the lower 
the feedback is, the quicker the stellar mass is assembled. 
In particular, throughout the first Gyr, FO builds a stellar 
component that is several times more massive than that of 
F80 (consistent with the results of S06), while the differ¬ 
ence between the runs with intermediate feedback and F80 
is less pronounced. This discrepancy slowly fades away with 
time: after 10 Gyr of evolution the total stellar mass assem¬ 
bled is virtually the same in all systems. Our interpretation 
is that, at the beginning of the simulation, stellar feedback 
quickly removes a significant fraction of gas that was piling 
up at the centre of the system, quenching star formation in 
the first few hundreds of Myr. A key point is that, in the 
feedback energy range explored in the models, this material 
remains gravitationally bound to the system, thus it slowly 
re-accretes onto the disc and takes part in the process of star 
formation. After 10 Gyr, the total amount of stars formed 
in the simulations with feedback and in that without feed¬ 
back is similar, as the total amount of fuel available does not 
differ. As FO has used most of its star forming fuel at the 
beginning of its evolution, the amount of neutral hydrogen 
left at 10 Gyr is 3.8 x 10® M©, ~ 13% of that available in 
the disc of F80 at the same time. 

It must be pointed out that the parameter F^sn, on 
which our analysis is focussed, is only one of the many in¬ 
gredients that determine the global ‘effectiveness’ of stel¬ 
lar feedback in simulations. For instance, [Governato et alT] 
(2010) realized the importance of the density threshold of 


gas eligible for being star-forming in simulations of dwarf 
galaxies: the larger the threshold, the more efficient the feed¬ 
back is in preventing the formation of a prominent stellar 
bulge by producing strong outflows from the center that re¬ 
move low angular momentum material. To explore this sce¬ 
nario, we re-ran our simulation F40 by increasing the star 


formation density threshold from our ‘fiducial’ value of 0.1 
cm“® (S06) to 100 cm“®. We found that the final (t = 10 
Gyr) stellar mass of this new system is similar to that re¬ 
ported in Table but the central stellar surface density is 
about 40% smaller than that of our fiducial run F40, indicat¬ 
ing that now feedback is slightly more efficient in preventing 
a mass stockpiling at the system centre. The total HI mass, 
also, is slightly larger (3.3 x 10^ Mq). The change in kine¬ 
matics is however negligible, as the peak rotational velocity 
of this new system is only 10kms“^ lower than that of our 
hducial run. A detailed analysis of the impact of the den¬ 
sity threshold on our systems is beyond the purpose of our 
study, but these results may indicate that massive systems 
like those studied in this paper are relatively less affected by 
large variation of this quantity. 

Finally we would like to stress that, even though our 
runs adopt a relatively high F^sn, in state-of-the art cos¬ 
mological simulations even larger values are required or ad¬ 
ditional feedback mechanisms (like radiative feedback from 
young stars or AGN feedback) are included in order to match 
the observations ( Hopkins et al.|[2014 Schaye et al.||2015 ). 
Our system with highest feedback shows a layer of extra- 
planar HI which is less massive and extended than that ob¬ 
served in NGG 891, a galaxy with a similar star formation 


rate and comparable total Hi and stellar masses (Ooster- 
loo et al!]|2007 Fraternali et al.pOll ). In future studies, it 
will be interesting to explore values of F^sn larger than those 
considered in this work. 


5.2 The origin of the rotational lag 


In Sections |3.1| and |4.2| we discussed the morphology and 
the kinematics of the extra-planar cold gas of the model 
with highest feedback, F80. To our knowledge, this is the 
hrst time that a global hydrodynamical simulation of a disc 
galaxy shows the presence of a thick Hi disc produced by 
stellar feedback with kinematics in agreement with the ob¬ 
servations. Melioli et al. (2009) used adaptive mesh rehne- 


ment hydrodynamical simulations to follow the dynamical 
evolution of a galactic fountain how powered by ~ 100 OB 
associations in a Milky Way-like galaxy. Both the supernova 
rate and energy input used in their work are comparable to 
those adopted in model F80. Melioli et al. (2009) found that 


gas from the ISM gets lifted up to a few kpc above the disc 
and falls back approximately to the same radius, in agree¬ 
ment with our results. However, they could not reproduce 
the lag in rotational velocity that is observed in the HI halo 
of real galaxies, and suggested that the galactic fountain 
how interacts with extragalactic material inhowing onto the 
galaxy with low angular momentum to produce these lagging 
kinematics. The same conclusions were previously found by 


Fraternal! & Binney (2008) (hereafter FB08) using a dy¬ 
namical model of the galactic fountain. In our simulations, 
a natural source of gas inhow is the gas corona. Following 
Peek et al. (2008), the net (inhow-outhow) gas accretion in 
a spherical shell at a distance r from the system centre can 
be determined in the models as 


dM(r) 

dt 


n{r) 


= E 


Ar 


(5) 
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F80, gas accretion 



Figure 11. Net gas flow (defined in eq.[^ as a function of the 
distance from the system centre derived at four different times 
in our model F80. Positive values indicate accretion. The solid 
lines represent pristine gas only, the dashed line (shown only at 
t = 10Gyr) includes polluted material as well. 


by evaluating eq.® ir 
at a distance 10°'^ 


where rrii is the mass of the z-th gas particle, Vr^i is the radial 
(i.e. directed towards the centre) component of its velocity, 
Ar is the thickness of the shell and the sum is extended 
to all n{r) particles inside that shell. Fig. 11 shows the net 
(inflow-outflow) gas flow as a function of the distance from 
the system centre in F80, derived at four different times 
in a series of shells, each one located 
kpc (0 ^ /c ^ 100) from the centre. At 
t = 10 Gyr, pristine gas (solid line) between 30 and 300 kpc is 
globally inflowing at a rate of about 2M0yr“^. At smaller 
distances, the inflow of pristine material appears to vanish 
as particles get contaminated by supernova feedback from 
the disc. However, when the total (polluted + pristine) gas 
is considered, the inflow appears to proceed down to the 
system centre (dashed line). In addition to the gas accreted 
from the corona, the discs of our systems are re-fuelled also 
by stellar winds, which return to the ISM about 25% of the 
gas converted into stars. Therefore, we estimate that the rate 
at which new gas is reverted to the disc of F80 at t = 10 Gyr is 
about 3 M© yr“^, which comes relatively close to the SFR of 
this system (4.2 M© yr“^). We also point out that in all the 
models (including FO) the mass flow profile is very similar, 
suggesting that stellar feedback does not have a big impact 
on the way the models are accreting gas from the GGM. 
Fig.[^ shows also that the gas accretion rate was higher 
at early times, peaking around 3.5M©yr“^ at t = 5 Gyr 
and 6M©yr“^ at t = 3 Gyr. Interestingly, this is in good 
agreement with the inflow rates derived from simulations of 
Milky Way-like galaxies in a full cosmological context (e.g. 
Brook et al.|2014 ). 


In model F80, it is not straightforward to determine the 
details of the interaction between the cold, metal-rich gas 
lifted from the disc and the hot coronal material inflowing 
to the centre, as pristine gas particles are absent in the re¬ 
gion where extra-planar cold gas exists. This is likely caused 
by stellar winds from evolved stars in the thick stellar disc, 
which inject mass and metals into the surrounding gas par¬ 
ticles (S06). In fact, pristine gas in F80 appears only above 
4 kpc from the disc, precisely where the stellar component 
fades out. This implies that, at the disc-corona interface. 


we can not anymore use the simple polluted/unpolluted cri- 
terium to determine to which component a gas particle be¬ 
longs. However, we found that the metallicity distribution 
of gas within a few kpc from the disc is bimodal, with two 
clear peaks located around Z = Z© and Z = 0.0IZ© and a 
separation occurring roughly at log(Z/Z©) = —1.2. For the 
purpose of investigating the disc-corona interplay, we inter¬ 
pret the gas at log(Z/Z©) > —1.2 as material ejected from 
the disc by stellar feedback and taking part in the galactic 
fountain gas cycle, and the gas at log(Z/Z©) < —1.2 as coro¬ 
nal gas which has never been part of the disc but is polluted 
by thick disc stars. Our interpretation is corroborated by the 
fact that gas particles at log(Z/Z©) < —1.2 are preferentially 
located outside the disc, while those at a higher metallicity 
are mainly located inside it. Fig.shows the rotation veloc¬ 
ity distribution for these two components in the region where 
extra-planar cold gas is present, 1 < \z\ < 4kpc, and com¬ 
pares it with that of the underlying Hi disc {\z\ <0.5kpc). 
Glearly, coronal gas rotates much slower than the material 
in the disc, whereas the galactic fountain gas distributes in 
between these two components. The kinematics of the coro¬ 
nal gas is not surprising since, given the high temperature, 
this material is partially pressure-supported and therefore 
must rotate at a velocity lower than that of the circular or¬ 
bit. Thus, a very likely interpretation is that metal-rich gas 
expelled from the disc interacts hydrodynamically with the 
slow-rotating and inflowing corona, loses angular momen¬ 
tum and produces a rotation velocity gradient in agreement 
with the Hi observations. This mechanism is qualitatively 
similar to that proposed by FB08. 

We stress that the loss of angular momentum experi¬ 
enced by fountain gas is likely to happen only at the disc- 
corona interface, where the specific angular momentum of 
the coronal gas is less than that of the disc material. How¬ 
ever, powerful feedback events can eject gas up to distances 
of several tens or hundreds of kpc, where this condition is 
reversed: this material can gain some angular momentum 
from the corona before being re-accret ed back to the disc 


(Brook et al. 


2012 


Ubler et al. 


2014). It is possible that 


such a mechanism occurs also in our runs since polluted gas 
is present up to 60 kpc from the system centres. 

The effects of the disc-corona interplay on the accre¬ 
tion of coronal gas have been studied in a series of papers 


(Marinacci et al. pon] [Marasco et al.|2012[ [Fraternali et al 


2013). In these works, a combination of dynamical models 
and high-resolution hydrodynamical simulations have been 
used to propose a scenario where the cooling of the lower 
corona is enhanced by stellar feedback from the disc via 
the galactic fountain mechanism. This process, known as 
‘supernova-driven gas accretion’, can hardly take place in 
our simulations as the resolution is not adequate to describe 
the turbulent mixing between these two gas layers. 


6 CONCLUSIONS 

We analysed the outcome of four N-body+SPH simulations 
of a Milky-Way like system produced by the radiative cool¬ 
ing of a rotating gaseous corona embedded in a dark matter 
halo. Each simulation uses the same initial conditions, but 
differs by the amount of energy per supernova, Esn, injected 
into the interstellar medium, which spans a factor of 32 from 
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F80, velocity distribution at the disc-corona interface 



Figure 12. Azimuthal velocity distribution of gas particles lo¬ 
cated at the disc-corona interface (1 < \z\ <4kpc) in our run F80, 
at t = 10 Gyr. The solid line shows the material ejected from the 
disc by stellar feedback, defined as having \og(Z/ZQ) >—1.2 (see 
text), while the dashed line shows the coronal gas. For compari¬ 
son, the dotted line shows the velocity distribution of cold gas in 
the disc. Each distribution is normalised to its peak value. 


the lowest feedback model (2.5 x 10^^ erg/SN) to the high¬ 
est feedback model (80 x 10^^ erg/SN). We studied how the 
morphology and the kinematics of these systems after 10 
Gyr of evolution compare to those of the Milky Way. In ad¬ 
dition, we performed a direct comparison with the Hi and 
absorption-line data of the Milky Way by simulating an ob¬ 
servation of the simulated systems from within the galaxy. 
Our main results are the following: 

• in the range of input energies considered, stellar feed¬ 
back has almost no impact on the star formation histories 
of the model galaxies; 

• stellar feedback significantly affects the vertical distri¬ 
bution of cold gas, as high feedback produces thicker discs 
with larger velocity dispersion; 

• the simulation with the highest feedback (F80) shows 
a prominent layer of extra-planar cold gas. The kinematics 
of this medium is consistent with that derived in real galax¬ 
ies like the Milky Way and NGC 891, although its column 
density is about one order of magnitude lower. This find¬ 
ing supports the idea that such features are produced by 
feedback activity from star-forming discs. 

• The location, kinematics and typical column density of 
the hot (5.3 < log(T) < 5.8) gas in F80 are consistent with 
those inferred in the Milky Way and external galaxies via 
OVI absorption line studies. Gas at lower temperatures is 
commonly observed in absorption around galaxies, but it is 
instead scarce in the GGM of our simulated systems. 

In addition, we found that the increase of Esn affects the 
stellar component in two ways: a) it systematically enhances 
the central stellar surface density, and consequently the con¬ 
centration of dark matter in the galaxy centres; b) it reduces 
the thickening and the velocity dispersion of the stellar disc. 
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